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ABSTRACT 

An approach was proposed and assessed for the high-fidelity modeling of 
progressive damage and failure in composite materials. It combines the Floating 
Node Method (FNM) and the Virtual Crack Closure Technique (VCCT) to 
represent multiple interacting failure mechanisms in a mesh-independent fashion. 
Delamination, matrix cracking, and migration were captured failure and migration 
criteria based on fracture mechanics. Quasi-static and fatigue loading were modeled 
within the same overall framework. The methodology proposed was illustrated by 
simulating the delamination migration test, showing good agreement with the 
available experimental data. 


INTRODUCTION 

Damage in composite materials generally occurs as a combination of different and 
interacting failure mechanisms, e.g. delamination and matrix cracking. Capturing 
these interactions accurately is essential to confidently model and predict 
progressive damage and failure. Several approaches have recently been proposed 
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that explicitly model different failure mechanisms and attempt to capture their 
interaction [1-3]. The present approach combines the Floating Node Method (FNM) 
[4] with the Virtual Crack Closure Technique (VCCT) [5,6] to explicitly account 
for different failure mechanisms and their interaction. Delamination, matrix 
cracking, and migration events are all modeled with the same FNM element using 
failure and migration criteria based on fracture mechanics. Preliminary results for 
quasi-static loading using this approach have been presented in [7]. The approach 
was validated using recent experimental results in which a setup capable of 
isolating a single complete migration event was developed [8], where migration is 
defined as the transition of a delamination from its current ply interface to a new 
interface via transverse ply cracking. In the present paper, some of the results 
obtained in [7] will be reviewed. In addition, the approach developed in [7] is 
combined with a fatigue algorithm and a novel migration criterion, and used to 
simulate delamination migration under fatigue loading conditions. 

FLOATING NODE METHOD (FNM) 

The FNM is a recently proposed numerical method, capable of representing 
multiple evolving discontinuities in solids [4], In the present section, an overview of 
the method is provided, further details and comparison with other existing methods 
are given in [4], 

Element formulation 

The static equilibrium of a body with volume ft under body forces with density f 
(acting on ft) and traction t acting on the boundary T n can be expressed in the weak 
form as: 


where u is the displacement vector; v is the test function; e is the strain tensor 
related to u through the differential operator relative to Cartesian coordinates L x as 
€ = £ x (u); and a is the stress tensor related to the strains through Hooke’s law as 
a = De, with D being the constitutive tensor. In the Floating Node Method, each 
real node “i” is characterized by its nodal coordinates x t and associated Degrees of 
Freedom (DoF) qj. In addition, an FNM element contains a suitable number of 
floating DoF without pre-defined associated nodal position vectors. These 
additional floating nodes are used to represent discontinuities. Their number varies 
with the number and type of discontinuities, weak or strong, modeled within each 
element. A strong discontinuity is defined as a jump in a field quantity (e.g. 
displacements), while a weak discontinuity is defined as a jump in the gradient of a 
field quantity (e.g. strains). Figure 1 shows an example of such an element with 
four nodes and four additional sets of floating DoF required to represent a strong 
discontinuity. 
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Figure 1. Overview of the Floating Node Method (FNM). 


ELEMENT FORMULATION WITH WEAK/STRONG DISCONTINUITY 


Once a discontinuity in the element is defined, the element is split in two or more 
partitions (depending on the discontinuity). Without loss of generality, a case in 
which the element is split in two partitions, U A and U B , is illustrated in fig. 1. For 
each partition, U A and U B , a vector of nodal coordinates, x^ A and x nfi , is defined. 
For the case in fig. 1 : 

X n A = [ x 5< x 6< x 3,X 4 ] (2) 


x n B = [xi,x 2 ,x 6 ,x 5 ] 
Each partition has a separate Jacobian: 

dx dN 

,A = d^ = df XnA 


(3) 

(4) 


dx dN 

,B = d5 = d5 x ^ (5) 

The displacements u A and u B , in partitions U A and U B , respectively, are 
interpolated separately from the DoF q A and q B associated with each partition: 

u A = Nq A and u A = Nq A (6) 

These DoF, q A and q B , did not have an associated position beforehand, and 
therefore were considered to be “floating”. As the discontinuity is defined, they are 
linked to the respective position vector. In the case of fig. 1, which represents a 
strong discontinuity, q A = [q5.q6.q3.q4] and q B = [q1.q2.q6.-q5.]- Note that 
there are four different sets of floating DoF: q 5 is different from q 5 , and q 6 is 
different from q 6 ,. If a weak discontinuity were to be modeled, only two sets of 




floating DoF would be included in the element, and the DoF with a prime would 
coincide with those without a prime. 

The strains then become: 


e A — ^x( u a) — £ 5 (N)J A 1 qA — B A q A 

(?) 

€ b = £x( u b) = A(N)J b 1( 1b = B B q B 
The stiffness matrices for partitions fi A and ft B are: 

(8) 

K a = BlDB A det(J A )dS 

(9) 

K b = £ B B DB B det(J B )dS 

(10) 

and the force vectors: 


Q a = l N T fdet(J A )dS + J N T f det(J A )dS 

(11) 

Q b = l N T fdet(J B )dS + J N T f det(J B )dS 

(12) 

Finally, the equations of equilibrium can be written as: 

K A q A = Q a and K B q B = Q B 

(13) 


ELEMENT TOPOLOGY AND ASSEMBLY 

To illustrate this approach, a floating node element capable of modeling two weak 
interfaces with an arbitrary crack between them was implemented, as shown in fig. 
2(a). The nodal position of the initial integration domain, and the real and floating 
DoF are indicated in the figure. Each floating DoF is associated with either an edge 
or the inner domain of the element. Once weak/strong discontinuities are detected, 
the floating DoF at the edges are used to determine the solution for the nodal 
positions created by the intersection of these discontinuities with the element edges 
(figs. 2(b) to 2(e)). The additional inner floating DoF are used to determine the 
solution for the nodal positions created by the intersection of multiple cracks within 
the element (fig. 2(f)). In figs 2(b) to (f), only the floating DoF used are 
represented. In each case, the floating DoF that are not used have no assigned 
position, only a topological relation to an edge or to the inner domain. During 
assembly, inner floating DoF can be removed from the analysis through static 
condensation. The floating DoF associated with the edges are assembled with the 
corresponding DoF of neighboring elements, and will therefore have a unique 
position in the global DoF vector. All floating DoF are then used as required to 
model weak/strong discontinuities as they evolve throughout the simulation. To 


accommodate certain crack geometries while enabling the integration of the domain 
of interest, the elements may need to be partitioned in triangles rather than 
rectangles only, figs. 2(d) and 2(f). Nevertheless, also in this case, the general 
procedure outlined in the previous section is used to integrate these elements, and 
determine their stiffness. 

The approach can be extended to model an arbitrary number of weak and strong 
discontinuities within an element, and therefore an arbitrary number of plies and 
interacting cracks, provided its topology, including the number of floating DoF, is 
adequately defined [4, 7]. 



(a) topology of the floating node element 
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Figure 2. Representation of the floating node element used in the present work, illustrating how the 
different weak/strong discontinuities scenarios are accommodated. 


VIRTUAL CRACK CLOSURE TECHNIQUE (VCCT) 


In the present work, the FNM method has been coupled with the Virtual Crack 
Closure Technique (VCCT) [5]. A comprehensive review of VCCT is provided in 
[6]. Traditionally, VCCT is used to obtain energy release rates in cases where the 
crack path is known beforehand and can be aligned with the elements’ edges, such 
as the case for delamination. In VCCT, the Mode I and II energy release rates are 
obtained from [6]: 


Gi 


1 

2 ct^ 



1/2 


(14) 








where F n and F t are the normal and tangential components of internal force vector F 
at the crack tip; and [q t ] are the normal and tangential components of the 
displacement jump [q] between the nodes immediately behind the crack tip; and 
and a 2 are the lengths of the crack in the elements behind and ahead of the crack 
tip, respectively (fig. 3). When combining VCCT with FNM, the forces and 
displacements needed to obtain the energy release rates, Eqs (14) and (15), are 
computed at the floating DoF as the crack develops (fig. 3). Additionally, floating 
DoF can also be added along a virtual crack plane up to a distance r, named the 
‘enrichment radius’, r enr . The associated floating nodes provide additional degrees 
of freedom along this path, which enable the accurate modeling of the deformations 
near the crack tip [4], At the end of the enrichment radius, the displacements of the 
floating DoF are interpolated from the displacements of the real DoF. In the present 
work, the enrichment radius was chosen to be equal to the largest in-plane 
dimension of the numerical model. 

A A A A A 



Figure 3. Virtual Crack Closure Technique and Floating Node Method for arbitrary crack 

propagation. 


DELAMINATION PROPAGATION 
Quasi-static loading 

In this paper, de lamination is modeled using YCCT to compute the energy release 
rates G, and G n at each delamination front position. These are then used in a failure 
criterion: 


G n 



(16) 


where G T = Gj + G n , and G c is the critical energy release rate given by [9]: 


G c 


— G lc + (G 


lie 



(17) 


For delamination growth according to this criterion, G c is assumed to equal the 
critical energy release rate of the interface, G l c nt . The delamination front is 
advanced by one element along the interface when G T = G l c nt . 


Fatigue loading 


Under fatigue loading, the delamination growth rate is determined using a Paris 
Law type relationship: 


where G Tmax = G Imax + G Umax . Since delamination is growing under mixed mode 
conditions, both A and n are assumed to vary with mode mixity. Several 
expressions have been proposed for characterizing the propagation rate under mixed 
mode conditions [10]. However, no general expression has yet found widespread 
acceptance. In the present work, a first order approximation is used to characterize 
the variation of A and n with mode mixity: 


where pairs A,, n h and A u , n u were obtained from fatigue testing under Mode I 
and Mode II conditions for the material being used, respectively [11, 12], Recent 
experiments indicate that a higher order approximation can possibly improve the 
representation of the variation of A and n with mode mixity [13]. Experimental 
results also confirm the monotonic behavior assumed in a first order approximation. 
Therefore, for simplicity, Equations 19 and 20 will be used in this study. 

The growth rate obtained from Equation 18 is then used to determine crack 
propagation using the fatigue algorithm implemented, which is discussed in a 
subsequent section. Note that Equation 18 does not account for the effects of 
loading ratio, and frequency. Therefore, its use is currently only recommended to 
predict the growth rate for R-ratios and frequencies that correspond to those used to 
generate the delamination fatigue characterization data. 

MATRIX CRACK PROPAGATION 
Quasi-static loading 

In the present study, matrix crack propagation is also modeled with VCCT. 
However, unlike delamination, the matrix crack is assumed to propagate following 
a Mode I path in the through-thickness direction, as supported by experimental 
evidence [8]. Therefore, the failure criterion is assumed to be given by: 


where G lc is assumed to be identical to the Mode I delamination fracture toughness, 
as demonstrated in [14], At each crack growth increment, the Mode I crack path 
orientation is approximately determined using a maximum tangential stress 
criterion. This criterion, typically written in terms of the Mode I and II stress 
intensity factors [15], can also be written in terms of energy release rates G, and G n , 



(18) 



(19) 


( 20 ) 


f(G T ) = ^ L - 1 = 0 


( 21 ) 


where the angle 6 that maximizes the tangential stress a ee is obtained by 
evaluating: 



( 22 ) 


and for the two solutions obtained choosing the one that maximizes a ee given by: 

{3 cos (£) + cos (y)] + ^ ^ 

/G^Sgn(F t ) |— 3 sin - 3 sin 

where F t is the tangential component of the internal force vector F at the crack tip. 
The tangential stress a ee relates to d^g by: 


where, 




a ee 



plane stress 
plane strain 


(24) 


(25) 


E is the Young’s modulus; and v the Poisson’s ratio. Once the angle is determined, 
the criterion given by Equation 21 is assessed. If the criterion is met, the crack 
advances by one element along the projected crack path. 


Fatigue 

Matrix cracks are assumed to follow a Mode I path in the through-thickness 
direction, also in fatigue. In this case, the growth rate is simply obtained by: 

^ = Ai(G max ) n ' (26) 

and used to determine crack propagation in the fatigue algorithm implemented, 
which is discussed in a subsequent section. The growth rate of a matrix crack under 
Mode I loading is assumed to be well approximated by the growth rate of a 
delamination growing under the same loading conditions. This approach represents 
an extrapolation of the observations made in [14] for quasi-static loading conditions 
to the fatigue regime, and requires further experimental validation. Similar to 
Equation 18, Equation 26 does not account for R-ratio or frequency effects, and 
therefore its use is currently only recommended for R-ratios and frequencies that 
correspond to those used to generate the delamination fatigue characterization data. 


MIGRATION 


Criterion for delamination to matrix crack migration 

STATIC 


In the present work, the criterion for delamination migration proposed in [7] will be 
used. In this criterion, all variables needed to determine whether migration occurs, 
can be directly obtained from the FNM results at each interface crack position. 
Delamination migration is observed to be preceded by the creation of local micro 
cracks, oriented as a function of the principal stresses ahead of the crack tip, as 
discussed in [16]. If energetically favorable, these micro cracks can accumulate, 
propagating through one of the bounding plies (matrix crack), rather than 
accumulating and propagating along (or near) the interface (delamination), and 
migration is obtained [8]. The criterion proposed in [7] uses the sign of the 
tangential component of the internal force vector F t , associated with Mode II shear 
displacement, to determine the direction of the micro cracks and therefore the 
material into which the interface crack would tend to migrate. Additionally, it 
assumes that migration can only be completed if it is energetically favorable. The 
criterion is given by M s \ { Gjf lt , Gf nt , F t } — > {—1,0,1} as: 



Sgn(F t ), 


Sgn(F t ), 


Sgn (F t ) < 0 


Sgn(F t ) > 0 


(27) 


where M s = 0 if the interface crack does not migrate, and M s = —1 or M s =1 
if the crack migrates into material A or B, respectively (see fig. 4). In the present 
form, Equation 27 assumes that the internal force vector is computed at the lower 
surface (Material B). If the internal force vector is computed at the upper surface 
(Material A), the signs of the inequalities in Equation 27 need to be reversed. 
Depending on the sign of F t , the ratio between the energy release rate determined at 


Qlnt Qlnt 

the interface, and the critical energy release rate of Material A or B, or 

Gn G r 


IS 


compared to the ratio between the energy release rate determined at the interface 

Qlnt Qlnt Qlnt 

and the critical energy release rate of the interface, If the ratio or is 

G c G c G c 

Qlnt 

greater than -f^, the necessary condition for migration is met, and M s = — 1 or 


M s = 1. Note that although this criterion is necessary, it is not sufficient for 
migration to occur unless the ratio being assessed is greater than one, i.e. the failure 
criterion for delamination or matrix cracking need to be met. Figure 4 illustrates 
how the migration criterion interplays with the failure criteria for delamination and 
matrix cracking, to determine whether delamination, migration or neither will 
occur. Figure 4 also shows the principal stress resultant for positive and negative 
applied shear, see fig. 4a and 4b, respectively. 
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Figure 4 Illustration of the combination between failure and migration criteria used to determine 
delamination migration, Equations 16 and 21, and Equation 27, respectively. 


FATIGUE 


The criterion used to predict migration under fatigue shares the same basic premises 
as the migration criterion used in the quasi-static case. The shear sign dictates the 
direction towards which the delamination will tend to migrate. To determine 
whether migration will occur, the fatigue crack growth rate in the material into 
which the delamination will tend to migrate is compared to the fatigue crack growth 
rate of the delamination. Similar to the static case, the criterion is given by 
M f : {GIF, G' nt , F t ) -» {-1,0,1} as: 
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Sgn(F t ), 


Sgn(F t ), 


Sgn(F t ) < 0 


Sgn(F t ) > 0 


(28) 


If Mf = 0 , the interface crack does not migrate, and if Mf = — 1 or Mf = 1 , the 
crack migrates into material A or B, respectively. The criterion is illustrated in 
Figure 5, for positive and negative shear sign applied. In the present work, the 
interfacial crack, or delamination, propagates between a 0° and a 90° ply. 
Assuming Material A is the 0° ply, and Material B is the 90° ply, the growth rate of 


a crack growing into a 0° ply is considered to be much lower than the growth rate 
for a crack to propagate at the interface, since fiber fracture is required for that 

process to occur. Thus, ( — ) is always much smaller than ( — ) . Therefore, even 

VdN/^ VdN/ /nt 


if the sign of F t favors migration into Material A (0° ply), according to Equation 28, 
migration will not occur. This criterion translates the hypothesis that, due to the 
difference in growth rates, the crack will advance along the interface before enough 
cycling has occurred for it to migrate towards Material A. In the present study, it is 
considered that the growth rate of a Mode I matrix crack can be characterized by 


Equation 26, while the delamination growth rate, 



, is assumed to be given 


by Equation 18. Since both are of the same order of magnitude, if the sign of F t 
favours migration into Material B (90° ply) depending on the relative magnitude of 

( — ) and ( — ) , migration can occur. 

VdN/^ \dN/ Int 



a) positive shear 



Mf = 0, interface crack 
propagation 


Mj-= 1 , migration 
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Int 


b) negative shear 

Figure 5 Illustration of the fatigue migration criterion, Equation 28. 


Matrix crack to delamination transition 


When a matrix crack reaches a weak interface, it is assumed to trigger 
delamination, both under static and fatigue loading conditions. This delamination is 
at first contained within the FNM element, as illustrated in fig. 2(f). In the 
following steps, the delamination will propagate, or migrate, as detailed in the 
previous sections. 




FATIGUE ALGORITHM 


The algorithm used to propagate fatigue cracks in the present work is illustrated in 
fig. 6, and is described next. 


For each crack tip, i calculate maximum Mode I 
and Mode II energy release rates. 


Using O’ determine the growth rate for each crack 
tip from the Paris Law, Eqs. 18, or 26. 


Using © determine the cycles/UV|* nc , needed to 
propagate each crack, i, by one element, Sa\ x el . 

From (§|, determine the minimum number of cycles 

required to propagate a crack over one 

element, 5 N \ ™ nc , and the correspondent crack tip n . 

Propagate the crack tip,n , determined in @ by 
one element. 

Update the cycle count: increase the total cycle 
count, N , by SN |™ nc ; accumulate SN\™ nc cycles to 
all cracks that did not propagate in this increment; 
initialize the number of accumulated cycles for the 
crack tip,n, that propagated in this increment. 



I | Finite element analysis 

| | Post-processing 

Figure 6 Fatigue algorithm implemented. 


In the first step, G Imax \ l and G IImax \ l are determined for each crack tip, i. In the 
present work, only the maximum energy release rate is used to obtain delamination 
and matrix crack growth rate, Equations 18, and 26. Often Paris Law equations, 
accounting for R-ratio and frequency effects, require the calculation of both G Imin 
and G IImin , or AG. In the context of the present algorithm, G Imin , G IImin could be 
obtained through an additional finite element analysis per increment. Having 
determined G Imax \ l and G IImax \\ the number of cycles, 8N \\ nc , that are needed to 
propagate each crack i by the length correspondent to the next element are obtained 
from: 


MIL 



(29) 


where Sa\ l le i is the length of the element ahead of the crack i, and SN \ l acc is the 
number of cycles accumulated by the crack tip before it grows to the next element. 
Subsequently, the minimum number of SN \ l inc and the correspondent crack tip are 
determined, SN \f nc and n, respectively. This determines the cycle increment in this 



iteration. The crack tip n is advanced by the correspondent 8 a |” ei . The accumulated 
cycle count of all crack tips (except n) as well as the total cycle count, are 
incremented by 8N |f nc . Finally, the accumulated cycle count for the crack tip n is 
set to zero, and a new increment starts. 

SIMULATION OF DELAMINATION MIGRATION SPECIMENS 
Delamination migration test configuration 

In reference [8], a test aimed at investigating delamination migration was proposed. 
The Delamination Migration (DM) test configuration is illustrated in fig. 7. The 
specimen has three key features that enable the controlled observation of 
delamination growth followed by migration to another ply interface. First, the 
specimen geometry is in the form of a beam with the intent of promoting uniform 
delamination growth and migration across the specimen width. Second, the 
specimen contains a Teflon insert (acting as an artificial de lamination) at an 
interface between a 0° ply (specimen span direction) and a stack of four 90° plies 
(specimen width direction) (fig. 7). This provides an opportunity for the 
delamination to migrate to another ply interface by kinking through the 90° ply 
stack. Third, the specimen can be loaded in a manner to cause delamination growth 
from the Teflon insert that eventually migrates to another ply interface. This 
sequence of fracture events is made possible by the way in which specimen loading 
affects shear stresses acting across the delamination front. 



Figure 7 Delamination Migration (DM) test configuration [8]. 

Finite Element Model 

The simulations presented in this work were all performed using the finite element 
solver Abaqus/Standard® 6.13 (Implicit). Plane-strain and a small-displacement 
formulation were used. The floating node element (fig. 2) was developed and 
implemented as a user defined element (UEL) in Abaqus/Standard®, and applied in 
the center region of the model, as shown in fig. 8. The specimen layup used in the 
experiments and in the model, is also given in fig. 8. In the layup description ‘T’ 
represents the Teflon insert. Each block of plies of the same orientation was 
modeled with a separate element CPE4 (through-thickness), except for the 
highlighted plies at the center of the model (fig. 8). These plies were modeled 


within a single FNM element, as illustrated in fig. 2. Rigid contact with friction is 
assumed between the bottom surface of the specimen and the baseplate, and the top 
surface of the specimen and the clamps. All the material properties used, including 
the friction coefficient, are given in Tables I, II and III. A displacement was applied 
at a single node, at different load offsets L, to simulate a displacement-controlled 
test, see fig. 8. All tests were performed under ambient laboratory conditions. 



Figure 8 Finite element model used. 

TABLE I. ELASTIC PROPERTIES, IM7-8552 [17] 


E ii 

^22 — ^33 

v 12 — v 13 

v 23 

g 12 — g 13 

^23 

161.0 (GPa) 

11.38 (GPa) 

0.32 

0.44 

5.17(GPa) 

3.98(GPa) 


TABLE II. FRICTION COEFFICIENT CARBON/EPOXY TO ALUMINUM [18] AND 
CRITICAL EN ERGY RELEASE RATES FOR DELAMINATION OF IM7-8552 [19]. 



d/ c 

G//c 

V 

0.23 

0.21 (kj/m 2 ) 

0.77(kj/m 2 ) 

2.1 


TABLE III. PARIS LAW COEFFICIENTS FOR MODE I AND MODE II OF IM7-8552 [11,12], 


Ai 

n, 

An 

n n 

8.757E-6 

6.71 

6.84E-7 

5.45 


BOUNDARY CONDITIONS 

As illustrated in fig. 7, the DM specimen is clamped at both ends. Each clamp was 
tightened using three screws at a fixed applied torque. The clamping conditions 
were simulated in a first step, by applying the clamping force, estimated to be 1700 
N, via two reference nodes coupled to the top surface of each clamp through a 
displacement constraint. Finally, and as was mentioned previously, rigid contact 
with friction was used to model the interaction between clamps, specimen and 
baseplate. To assess the adequacy of the boundary conditions applied, the deflection 
of a DM specimen was measured using Digital Image Correlation (DIC) and 
compared to the numerical results. Figure 9 compares the DIC results with 
deflections obtained from three different numerical models. The curve labeled 
‘FIXED’, was obtained using a numerical model with nominal dimensions, and 


where the clamping conditions were idealized by fixing the displacement at the top 
and bottom surfaces of the specimen in the clamped region. It is evident that this 
model under-predicts the deflection. Furthermore, the maximum load is also not 
aligned with the maximum deflection observed experimentally. Measuring the exact 
load-application point and crack-tip position of the specimen tested, it was observed 
that both deviated from the nominal values. Including this correction in the model 
the curve ‘FIXED-CORR’ was obtained. The maximum deflection in this curve is 
now aligned with the observed experimentally, but the magnitude is still 
underestimated. Finally, the third curve, ‘CLAMP’, is obtained by applying the 
clamping force to the model in an initial step (including the corrected crack length 
and load-application point), and modeling the contact between clamps, specimen 
and baseplate, as described above. It is evident that these boundary conditions lead 
to the best agreement with the experimental data. However, the deflection seems 
slightly overestimated, near the left-hand-side clamp, which could be caused by an 
underestimation of the local clamping force, or friction coefficient. 



x, mm 

Figure 9 Comparison between a specimen deflection measured using DIC and the simulations. 

Results 

QUASI-STATIC LOADING 

Overall, the framework proposed is capable of predicting and simulating 
delamination migration in cross-ply laminates under quasi-static conditions. Fig. 10 
compares the computed load-displacement curves obtained, for two values of L, 
with experiments illustrating the accuracy of the approach. Overall, the maximum 
load is well predicted for both cases. For L = a 0 , the simulations predict unstable 
delamination after the peak load (fig. 10(a)). This unstable delamination stops just 
before migration. After the migration event, delamination arrests, and upon further 
loading it continues propagating stably along the next 0°/90° interface. A similar 
sequence of events was observed experimentally [7]. However, in the experiments, 
after the first unstable event, further loading was necessary before migration 
occurred. For the case L = 1.3a 0 , the simulations predict that the first peak in load 
is followed by a region of stable delamination growth. A region of unstable 


delamination, corresponding to the load-drop follows as shown in fig. 10(b). The 
migration event is predicted to occur in this region. After the migration, unstable 
delamination continues propagating along the next 0°/90° interface. Loading the 
specimen further eventually leads to a transition from unstable to stable 
delamination propagation. This sequence of events is identical to what was reported 
experimentally [8]. The only difference observed is that in the experiments, after 
the peak load, delamination propagates through a series of unstable events, rather 
than stably, as predicted by the simulations. Further simulations for other load- 
offsets, namely L = l.la 0 and L = 1.2 a 0 were performed and compared to the 
experimental results, showing a similar agreement. 




Figure 10 Load displacement curves for (a) load-offset L=a 0 and (b) L=1.3a 0 

FATIGUE LOADING 

Constant amplitude fatigue conditions were simulated. A maximum displacement 
corresponding to G Tmax = 0.8G C was applied, where G c was determined from the 
quasi-static simulations. This G c corresponds to the critical energy release rate to 
initiate delamination for each load-offset L. The frequency of the fatigue loading 
was assumed to be 5Hz, and the ratio of minimum to maximum displacement, R- 

ratio, ^ 2 ^ mtn = o.l. Both R-ratio and frequency are the same used to obtain the 

fatigue characterization data summarized in Table III. 

Overall the predicted damage morphology was similar in both quasi-static loading 
and fatigue, only fatigue results are shown here for brevity. Figure 11(a) provides 
an illustration of a typical crack path predicted by the simulations. Figure 11(b) 
shows the distance from the load-application point to migration, A K , as a function of 
the normalized load offset, L/a 0 . The simulations predict a decreasing trend in A K 
with L/a 0 , which is similar to what was observed in the static case [8]. 

Figure 12(a) shows two curves corresponding to the number of cycles to migration 
onset and to 6a = 3 mm past the completed migration as a function of the load- 
offset, labeled ‘MIGRATION - INI’, and ‘MIGRATION + 8a\ respectively. 
Comparing ‘MIGRATION - INI’ to ‘MIGRATION + da’, it is clear that they only 
differ noticeably for L = a 0 . Overall, simulations predict a significant increase in 
cycle count to migration as the load-offset increases. This trend can be further 
examined by analyzing the crack growth as a function of number of cycles. Figure 
12(b) presents delamination growth up to migration for two different load-offsets 


L = a 0 and L = 1.3a 0 . Comparing the results obtained for the two offsets, it is 
possible to see that for L = 1.3a 0 a significant decrease in crack growth rate can be 
observed, between N=50000 and N=1 50000, before it increases again leading to 
further delamination and migration. 


U- 



6 

A k , 

mm 

3 


0 4 T T 1 

1 1.1 1.2 1.3 

Normalized load offset, L/a 0 

Figure 1 1(a) typical crack path, (b) distance from the migration location to the load-application 

point, A k , as a function of load-offset L. 
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Figure 12 (a) number of cycles to migration onset, and to migration onset + 6a = 3 mm, and (b) 
delamination growth prior to migration as a function of number of cycles. 


DISCUSSION 

The detailed analysis of the specimen deflection using DIC has helped in assessing 
the adequacy of the boundary conditions applied. Nevertheless, the clamping force 
and friction coefficient are likely to vary slightly from specimen to specimen. 


Therefore, further sensitivity studies are recommended to characterize fully the 
effect of the boundary conditions on the overall response of the specimens. 

Overall, the sequence of events, failure morphology and load-displacement curves 
are in good agreement with the experimental results. However, for the case L = a 0 , 
experimental results show an increase in load at the end of the first unstable event, 
before migration, not predicted by the simulations. Also, for L = 1.3a 0 , the 
simulations predict stable delamination, prior to the final load-drop, while 
experiments show that the delamination propagates through a series of unstable 
events. One of the main differences between simulations and experiments, is the 
assumption that the delamination propagates at the interface between the 0° and 90° 
ply. However, in the experiments, it is often observed that delamination delves 
slightly into 0° ply propagating just within it, and leads to the development of 
bridging fibers. This mechanisms is thought to be caused by the applied shear force 
that tends to drive the delamination into the 0° ply in the first part of the test (for the 
load offsets tested). Delamination delving, and the associated fiber bridging, leads 
to local variations in G c not accounted for by the model, and likely to lead to some 
of the discrepancies observed. 

The predicted damage morphology under fatigue loading conditions is similar to the 
one obtained under quasi-static loading, including migration location [8]. Overall, 
simulations predict an increase in cycles to migration as a function of the load- 
offset. For the L = 1.3a 0 case the delamination growth rate decreases substantially 
before it increases again and migration occurs. This decrease indicates that for 
lower displacements applied, a fatigue threshold may be reached, before migration 
is obtained. Further, also the delamination delving mechanism observed in the static 
tests might occur in fatigue, since it is a function primarily of the shear sign applied. 
The delving may alter (decrease) the growth rates obtained in the characterization 
tests. Hence, threshold values might be reached even for high-applied loads, and 
overall the results obtained with the present approach, may be overly conservative, 
i.e. more crack growth is predicted for the same number of applied cycles. In the 
present work, the fatigue loading simulated was chosen to be consistent with the 
characterization data available, i.e. same R-ratio, and frequency. A more general 
fatigue loading scenario would require either further characterization data, or the 
use of a modified Paris Law that can potentially account for these effects [20]. 
Furthermore, the first order approximation used to determine the variation of 
delamination growth rate with mode-mixity, Equations 19 and 20, can potentially 
be improved given the recent experimental data obtained [13]. Finally, although the 
fatigue algorithm and overall approach has been demonstrated, experimental 
validation is needed. 


CONCLUSIONS 

The FNM (Floating Node Method) was combined with YCCT (Virtual Crack 
Closure Technique) and used to model delamination migration in cross-ply 
laminates under quasi-static and fatigue loading conditions. Delamination, matrix 
cracking, and migration, are all modeled with the same FNM element using failure 
and migration criteria based on fracture mechanics. As a result, the present 
approach only uses material properties, the majority of which can be obtained using 


well established standard characterization tests. For quasi-static loading, the results 
obtained compared favorably with experimental results available. Regarding 
fatigue, the algorithm and migration criterion have been demonstrated. However, 
experimental validation is still required. 
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